//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark all sample
reg $baseline_child
gen sample0 = e(sample)

// mark baseline sample
reg $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen sample1 = e(sample)

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'																									   if sample1 == 1
	gen  `var'_std = (`var' - r(mean))/r(sd)                                                                       if sample1 == 1
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std) if sample1 == 1, method(adf);
predict parenting_age1 if sample1 == 1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std) if sample1 == 1, 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3 if sample1 == 1, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3) if sample1 == 1

global inputs_std   cum_avg_daycare_36m_sum parenting_ages13
keep   ihdp site tg bwg sample* $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs ppvtstd wasifsiq

// merge in ages 5 and 8
merge 1:1 ihdp using iqchildhoodihdp
keep  if _merge == 3
drop _merge

// impute iq
replace kidppvt5 = wippsif5 if kidppvt5 ==.
replace wippsif5 = kidppvt5 if wippsif5 ==.

replace ppvtstd8 = fsiq8    if ppvtstd8 ==.
replace fsiq8    = ppvtstd8 if fsiq8    ==.

replace ppvtstd  = wasifsiq if ppvtstd  ==.
replace wasifsiq = ppvtstd  if wasifsiq ==.

reg fsiq8 ppvtstd8 kidppvt5 wippsif5 wasifsiq ppvtstd
gen sample2 = e(sample) 

// outcomes
merge 1:1 ihdp using adultihdp
keep if _merge == 3
drop    _merge

// construct "positive" outcomes
replace idle18y         = 1 - idle18y
replace sch_eversped18y = 1 - sch_eversped18y
replace sch_math18y     = 1 - sch_math18y
replace sch_read18y     = 1 - sch_read18y
replace sch_thrp18y     = 1 - sch_thrp18y
replace teen18y         = 1 - teen18y

gen     smoke18 = .
replace smoke18 = 1 if cigs18y >  0  & cigs18y != .
replace smoke18 = 0 if cigs18y == 0 & cigs18y != .
replace smoke18 = 1 - smoke18

gen     absence18 = .
replace absence18 = 0 if sch_dabs18y < 5  & sch_dabs18y != .
replace absence18 = 1 if sch_dabs18y >= 5 & sch_dabs18y != .
replace absence18 = 1 - absence18 

// marking missing
global ed18  sch_eversped18y sch_math18y sch_read18y sch_test18y
global beh18 smoke18 idle18y sch_thrp18y teen18y
egen missing     = rowmiss($ed18 $beh18)

// average
egen avg_outcome_educ18 = rowmean($ed18)        if missing == 0
egen avg_outcome_beh18  = rowmean($beh18)       if missing == 0
egen avg_outcome18      = rowmean(avg_outcome_educ18 avg_outcome_beh18) if missing == 0

// define sample
reg     avg_outcome18 if sample2 == 1
replace sample2 = 0   if e(sample) == 0
replace sample2 = 0   if   sample1 == 0

// lists of outcomes
global iq5    kidppvt5 wippsif5

// describe by sample
global sample_all          if twin != .
global sample_singletons   if twin == 0
global sample_twins        if twin == 1

// describe
cd $output
foreach group in all singletons twins {
	foreach var of varlist kidppvt5 wippsif5 {
		bootstrap, strata(bwg site tg) reps($bootstraps) saving(b`var'_`group', replace) : reg `var' tg ${sample_`group'} & sample2 == 1
		matrix b`var'_g`group' = e(b)[1,1]
	}
	matrix db_g`group' = bkidppvt5_g`group' - bwippsif5_g`group'
	
	foreach var of varlist kidppvt5 wippsif5 {
		preserve
		use  b`var'_`group', clear
		rename _b_tg b`var'_`group'
		gen n = _n
		save b`var'_`group', replace
		restore
	}
	preserve
	use bkidppvt5_`group', replace
	merge 1:1 n using  bwippsif5_`group'
	gen    db_g`group'B = bkidppvt5_`group' - bwippsif5_`group'
	summ   db_g`group'B
	matrix dbse_g`group' = r(sd)
	restore
	
	matrix  dbt_g`group' = abs(db_g`group'[1,1]/dbse_g`group'[1,1])
	matrix  dbp_g`group' = 2*(1 - normal(dbt_g`group'[1,1]))	
	matrix `group' = [db_g`group',dbse_g`group',dbp_g`group'] 
}
mat list singletons
mat list twins
mat list all 
